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The objective of this study is to conduct a unified computational analysis for computing 
design parameters such as axial thrust, convective and radiative wall heat fluxes for 
regeneratively cooled liquid rocket engine nozzles, so as to develop a computational strategy for 
computing those parameters through parametric investigations. The computational methodology 
is based on a multidimensional, finite-volume, turbulent, chemically reacting, radiating, 
unstructured-grid, and pressure-based formulation, with grid refinement capabilities. Systematic 
parametric studies on effects of wall boundary conditions, combustion chemistry, radiation 
coupling, computational cell shape, and grid refinement were performed and assessed. Under the 
computational framework of this study, it is found that the computed axial thrust performance, 
flow features, and wall heat fluxes compared well with those of available data and calculations, 
using a strategy of structured-grid dominated mesh, finite-rate chemistry, and cooled wall 
boundary condition. 

Nomenclature 

C], C 2 , C s , C turbulence modeling constants, 1.15, 1.9, 0.25, and 0.09. 

C p = heat capacity 

D = diffusivity 

H = total enthalpy 


* Staff Consultant, Applied Fluid Dynamics Analysis Group, TD64, Senior Member AIAA. 
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h = static enthalpy 

I = radiative intensity 

K = thermal conductivity 

k - turbulent kinetic energy 

P = pressure 

Q = heat flux 

R = recovery factor 

r = location coordinate 

T = temperature 

T* = law-of-the-wall temperature 

t = time, s 

u, v, w = mean velocities in three directions 

«r = wall friction velocity 

x = Cartesian coordinates 

a = species concentration 

e = turbulent kinetic energy dissipation rate 

6 = energy dissipation contribution 

k = absorption coefficient 

H = viscosity 

Ht = turbulent eddy viscosity (=pC |1 k 2 /e) 

II = turbulent kinetic energy production 

p = density 

a = turbulence modeling constants 
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x = shear stress 

Q = direction vector. £T denotes the leaving radiative intensity direction 

0 ) = chemical species production rate 

Subscripts 
b = black body 

c = convective 

cl = centerline 

p = off-wall (wall-function) point 

r = radiative 

t = turbulent flow 

w = wall 

0 = reference 


I. Introduction 

The two major factors in rocket engine design, performance and integrity (convective heat 
transfer), are often conducted separately. As a result, the final design based on performance may 
have to be altered due to convective heating considerations, and vice versa, resulting in delays 
and compromises. Recently, radiative heating has been generating concerns due to renewed 
interest in hydrocarbon engines. Those combined motivated this study to perform and 
demonstrate a unified analysis for the computation of those design parameters in a simultaneous 
fashion. Systematic parametric studies on effects of wall boundary conditions, combustion 
chemistry, radiation coupling, computational cell shape, and grid refinement were performed and 
assessed, in order to come up with a strategy for efficient and realistic computations. 
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An axisynunetric nozzle axial force analysis, 1 and a conjugate convective heat transfer 
analysis 2 for the Block I Space Shuttle Main Engine (SSME) thruster were reported in the 
1990’s, using a structured-grid, multi-zone computational fluid dynamics (CFD) solver FDNS. 
As the requirements for parallel computing efficiency and faster grid generation arise, an 
unstructured-grid CFD methodology UNIC was developed recently through several activities, 
namely the launch vehicle base-heating, 3 Laser propulsion, 4 and stage-separation. 5 This 
unstructured-grid CFD methodology is refined in this study to conduct a series of unified axial 
force, convective and radiative heat transfer analyses, simulating SSME hot-firing at sea level. 
Both axisymmetric and three-dimensional analyses were performed and a final computational 
strategy reported. 


11. Computational Methodology 

The time-varying transport equations of continuity, species continuity, momentum, global 
energy (total enthalpy), turbulent kinetic energy, and turbulent kinetic energy dissipation can be 
written as: 




dpCCi 


+ _d 

dt dx 


-{fiUjOCj )= 


dpuj 

dt 


j 

a 




pD+ 


Mt 




'ctj 




+<y. 


dxj 


Wj“i)=- 


dx t dxj 


i£ h )^ q *((JL + JL 

dt dx j J dt dx j cr H 


w 


ff 


d x j 


(m+m,)- 




\Cp a H 


y{v 2 i2j 


We 


4 



d pk { d 

dt dxj 

dp£ d 
dt dxj 


■{oujk)= 

(puje)= 


dxj 


dxj 


a k) 


dk 


fl + 






'ej 


dxj 

de 


dxj 


+/K rw) 

+pj(c l n-c 2 €+c 3 n 2 /e) 


where the energy dissipation contribution 0 can be expressed as: 
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and the shear stress tjj can be expressed as: 
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A predictor and corrector solution algorithm was employed to provide coupling of the fluid 
governing equations. A second-order central-difference scheme was employed to discretize the 
diffusion fluxes and source terms of the governing equations. For the convective terms, a 
second-order upwind total variation diminishing difference scheme was used in this effort. To 
enhance the temporal accuracy, a second-order backward difference scheme was employed to 
discretize the temporal terms. A point-implicit (operator splitting) method was employed to solve 
the chemistry system. 


An extended k-e turbulence model 6 was used to describe the turbulence. A 7-species, 9- 
reaction detailed mechanism 7 was used to describe the finite-rate, hydrogen/oxygen (H 2 /0 2 ) 
afterburning chemical kinetics. The seven species are H 2 , 0 2 , H 2 0, 0, H, OH, and N 2 . 

A modified wall function approach was employed to provide wall boundary layer solutions 
that are less sensitive to the near-wall grid spacing. Consequently, the model has combined the 
advantages of both the integrated-to-the-wall approach and the conventional law-of-the-wall 
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approach by incorporating a complete velocity profile and a universal temperature profile 7 . This 
approach is especially useful in three-dimensional (3-D) applications. 

The convective heat transfer follows the modified Newtonian law 

The radiative heat transfer is analyzed by solving the radiative transfer equation 

(ft • V)J(r,£2) = -Kl(r,Q.) + tcl b (r) 

Discrete ordinate method was used to solve the radiative transfer equation and H 2 O is the major 
radiating medium. The spectral-line based weighted-sum-of-gray-gases model 8 was used to 
calculate the total emissivity and absorptivity of the radiating medium. Details of the numerical 
algorithm can be found in Refs 3-5. The radiative heat flux is given by the integration of the 
wall leaving radiative intensities 

Q r = I 7(r,a")|nT2 _ |dQ" 

n.Q _ <0 

III. Computational Grid Generation 

The flowfields of four axi symmetric and four 3-D grids were computed during the course of 
this study. The results from two representative axi symmetric and two 3-D grids are reported for 
conciseness. These grids are hybrid grids and can be classified into two groups: the structure- 
grid dominated and the unstructured-grid dominated. Figure 1 shows the layout of an 
axisymmetric, unstructured-grid dominated hybrid grid axi. It has four layers of structured 
(quadrilateral) grid surrounding the solid walls, while the rest of the domain filled with 
unstructured (triangular) cells. These structured-grid layers are used to ensure proper wall 
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boundary layer development. The layout of an axisymmetric, structure-cell dominated hybrid 
grid ax6 is shown in Fig. 2. The interior of the thruster and plume region is filled with 
quadrilateral cells, while the rest of the domain is filled with triangular cells. The structured-grid 
layers used in grid ax 1 are also embedded in grid ax6, and in 3-D grids 3d6 and 3d9 (Figs 3-4), 
such that the boundary layer development for all grids is similar. Figure 3 shows the layout of 
the hybrid 3-D grid 3d6. It was constructed by rotating grid ax6 72 times for 360 degrees. 
Figure 4 shows the layout of the hybrid grid 3d9. These computational grids were generated 
using the software package GRIDGEN. 9 Table 1 shows the total number of points and cells in 
these four grids. The structured cells in grids 3d6 and 3d9 are hexahedral elements. The 
unstructured cells in grid 3d9 are tetrahedral elements, while the unstructured cells in grid 3d6 
are prismatic elements. Note that all four grids were computed as a single zone, thus avoiding 
the interface complexities commonly seen in multi-zonal grids. 


grid 

# points 

# cells 

# structured cells 

# unstructured cells 



30,578 

2,016 

28,562 


17,391 

17,710 

15,300 

2,410 


1,286,934 

1,275,120 

1,101,600 

173,520 


418,165 

1,732,081 

227,984 

1,504,097 


IV. Boundary Conditions and Run Matrix 

Fixed total condition was used for the inlet of the thruster and the outer boundary. A total 
pressure of 1 atm was specified for the outer boundary in order to simulate the nozzle hot-firing 
at sea level. No-slip boundary condition was specified for the thruster walls. Symmetry 
condition was applied to the centerline for axisymmetric cases. The CEC program 10 was used to 
obtain the chamber equilibrium species composition for use at the thruster inlet. 
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The run matrix is shown in Table 2. These cases were built up systematically in order to 
understand the grid effects such as cell shape and grid refinement, and the modeling effects such 
as chemistry, wall boundary condition, and radiation. For the convenience of presentation, 
abbreviated letters are used to represent different cases in the run matrix. For example, case fz 
represents parametric conditions of frozen chemistry and adiabatic wall, while case frcgr uses 
parametric conditions of finite-rate chemistry, cooled wall, with grid refinement and radiation 
coupling. Due to the limitation of current resources, grid refinement was not performed for the 
3-D cases. 


Table 2. F 

tun matrix 

case 

chemistry 

wall 

grid 

refinement 

radiation 

fz 

frozen 

adiabatic 

off 

off 

eg 

equilibrium 

adiabatic 

off 

off 

ff 

finite-rate 

adiabatic 

off 

off 

fro 

finite-rate 

cooled 

off 

off 

frcr 

finite-rate 

cooled 

off 

on 

frcg 

finite-rate 

cooled 

on 

off 

frcgr 

finite-rate 

cooled 

on 

on 


V. Results and Discussion 

The computations were performed on a cluster machine using four processors for each 
ax i symmetric case and thirty-two processors for each 3-D case. A global time step of 1 ps was 
used. Figure 5 shows a comparison of the computed adiabatic wall temperatures for grid ax6. 
Similar, if not identical wall temperature profiles were obtained for grid axl and are not shown. 
The computed wall temperature for the frozen chemistry case is constant, indicating the 
conservation laws are satisfied. Those for the equilibrium and finite-rate chemistry cases increase 
first after the throat, due to the recombination of chemical species to become H 2 0. The 
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temperatures for those two cases then decrease as H 2 O dissociates. Of interest is the temperature 
for the equilibrium case, it drops continuously until it closes to that of the frozen flow, near the 
nozzle exit. This is expected since the stagnation temperature is very close to the chamber 
temperature with which the frozen composition was determined with an equilibrium solution. 
This also implies the equilibrium chemistry probably dissociates the H 2 0 at too fast a rate inside 
the nozzle. A specified cooled wall temperature profile 1 is also shown in Fig. 5 which was 
determined through a separate conjugate heat transfer calculation; this temperature profile is used 
later as a boundary condition to consider the effect of heat loss to the regenerative coolant 
channels. 

Figure 6 shows the computed Mach number contours for cases fz and frcgr for grids axl and 
ax6, respectively; those for other cases are similar to those of case frcgr and are not shown. 
These figures show the captured nozzle flow features (nozzle shock, lip shock, triple point, Mach 
disc, shock reflection, and shear layer/shock interaction). In general, all cases capture the flow 
features reasonably well, except the frozen flow case in which a curved Mach disk was obtained. 
It can also be seen that the nozzle shocks appear to be sharper in the contours of the structured- 
cell dominated grid ax6 than those of the unstructured-element dominated grid axl. The sharpest 
nozzle flow features are captured with grid refinement on grid ax6, while the added radiation 
changes the flow features only slightly. 

The significance of a curved disk is that a large flow recirculation appears behind the curved 
disk. The occurrence of the curved disk may be attributed to the difference in thermodynamics 
between the frozen flow and chemically reacting flows. As shown in the centerline H 2 0 mass 
fraction, specific heat ratio and Mach number profiles in Fig. 7, the inability of recombining the 
species of the frozen flow results in much higher specific heat ratios than those of the reacting 
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flows, which in turn produces high Mach numbers along the centerline. The higher shock 
strength leads to higher total pressure loss across the shock, which causes the shock center to 
retreat and consequently an overall curved disk. On the other hand, the curves of equilibrium 
chemistry closely follow those of finite-rate chemistry. This is because the centerline 
temperatures drop continuously and are much lower than the chamber temperature (Fig. 8), 
hence the dissociation process occurring on the adiabatic wall is frozen on the centerline. 
Furthermore, the curved disk phenomenon happens both in grid axl and ax6, hence it is cell- 
shape independent and thermodynamically induced. 

Figure 8 shows a comparison of the thruster centerline temperatures for grid ax6. The frozen 
chemistry gives the lowest bound while all other cases group together as an upper bound, while 
the result from Ref. 1 falls in between. As discussed above, the low frozen chemistry curve is 
caused by the thermodynamics. As for the result from Ref. 1, it is speculated that an older 
thermodynamics database was used then. A comparison of the thruster wall pressures is shown 
in Fig. 9. The computed results from all cases appear to group together and agree reasonably 
well with the test data. Figure 10 shows a comparison of thruster centerline pressures. All 
predictions agree reasonably well, except for the frozen flow case that deviates lower near the 
nozzle lip. 

Figure 1 1 shows the computed convective heat fluxes for grid ax6. As expected, the peak 
convective heat fluxes occur at the throat (x = 0) for all cases. The refined grid gives a slightly 
lower peak heat flux. Radiation does not affect the convective heat flux, because the maximum 
radiative heat flux is about two orders-of-magnitude lower than that of convection (Fig. 12). All 
predictions compare reasonably well with those of the three other methods. 2,11 Result from Fig. 
11 demonstrates that both the momentum and thermal wall boundary layers were captured 
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reasonably well with the current methodology. The difference in the initial heat fluxes is caused 
by the difference in ways of initiating the boundary layers among different methods and the 
significance of which is negligible in comparison to the peak heat flux. 

Figure 12 shows the computed radiative heat flux for grid ax6 while cooled wall, finite-rate 
chemistry, and grid refinement were used as operating conditions. As expected, high radiative 
heating occurs inside the combustion chamber within which the high temperature and high H 2 0 
concentration are prevalent. As the propulsive flow expands past the throat, the temperature 
drops, hence the low radiative heat flux. The peak radiative heat flux is about two orders of 
magnitude lower than that of the convective heat flux, which is reasonable for a hydrogen fueled 
engine. In current methodology, the injector faceplate is modeled as a black body. In order to 
compare the predicted radiation with that of a plume radiation code GASRAD, 12 which does not 
model the injector faceplate, another run was performed by setting the temperature of the injector 
faceplate to 300 deg. K, effectively turning off the black body radiation. The structured-grid 
solution from Ref. 1 was used as input for GASRAD radiation calculation, since GASRAD can 
not read unstructured-grid information. It can be seen that the result from turning off the 
blackbody radiation at the inlet, using a weighted-sum-of-gray-gases (WSGG) absorption model, 
compares reasonably well with that of GASRAD in which a narrow band (NB) absorption model 
was used. It should be pointed out that GASRAD reads in flow solution for a decoupled 
radiation solution, and current methodology solves the flow equations and radiative transfer 
equation simultaneously. The computed peak value is higher when the black body radiation is 
included at the inlet, as expected. It is also noted that GASRAD was developed for the prediction 
of plume radiation, hence it does not consider the re-radiation from the solid walls. In addition, 
it solves the line-of-sight equation and not the radiative transport equation. 
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Ta ble 4. Comparison of SSME thrust chamber specific impulses (ISP’s) for axisymmetric cases 



axl 

ax6 

fz 

438.7 

439.7 

eq 

455.6 

456.0 

fr _ a 

455.2 " 

455.6 

frc 

452.4 

452.6 

frcg 1 

452.9 

453.7 

frcgr ; 

452.5 

453.3 

Data 

453 

1.3 


Table 4 shows the comparison of computed SSME thrust chamber specific impulses, or the 
axial thrust performances for the axisymmetric cases. The frozen flow calculations give too low 
an axial force, even with the adiabatic wall assumption that assumes zero wall heat loss. This is 
again caused by inadequate heat capacity distributions forced by a fixed species composition. 
The reacting flow (with adiabatic wall) cases overpredict the data for about 2-3 s, with the 
equilibrium case giving the highest values. When the wall heat loss is considered (case frc), the 
axial force predictions become very close to the data. The quadrilateral cell dominated grid ax6 
appears to predict slightly better ISP’s than those of the triangular cell dominated grid axl. 
Within grid ax6, the grid refinement and radiation options (case frcgr) provide the best 
agreement. 

Figure 13 shows the computed temperature contours for grid 3d6, case frc. Similar to the 
Mach number contours, the temperature contours also show the captured nozzle flow physics 
such as the nozzle shock, lip shock, triple point, Mach disc, shock reflection, and shear 
layer/shock interaction. Two perpendicular planes are used to give the Mach disc a three- 
dimensional feel. The high temperature in the mixing layer indicates afterburning. 

Figure 14 shows a comparison of computed thruster centerline temperatures. The centerline 
temperature of grid 3d6 matches that of grid ax6 reasonable well, except inside the chamber 
where the temperature of grid 3d6 is slightly lower. Figure 15 compares the wall pressures. The 
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wall pressures of grid 3d6 and ax6 overlap and both compare reasonably well with the data. 
Figure 16 compares the centerline pressures. The centerline pressure of grid 3d6 coincides with 
that of grid ax6, until the nozzle lip where the pressure of grid ax6 is slightly higher. 

Figure 17 shows a comparison of convective wall heat fluxes. The computed 3-D heat fluxes 
agree reasonably well with those of other methods, while the predictions of grid 3d6 overlap with 
those of grid ax6. The radiation does not affect the convective heat fluxes of grid 3d6, again due 
to the relative low radiative heat fluxes inside a H 2 /O 2 engine. Figure 18 shows a comparison of 
the computed radiative wall heat fluxes. Similar to the result of the axisymmetric cases (Fig. 
12), the computed 3-D radiative fluxes using a weighted-sum-of-gray gases absorption model 
compares reasonably well with that of GASRAD using a narrow band absorption model, when 
the blackbody radiation at the inlet is turned off; And the predicted radiative heat flux is higher 
while the blackbody radiation at the inlet is turned on. 


Table 5. Comparison of SSME thrust 



3d6 

3d9 

fz 

439.6 

436.5 

_eg_ 

456.4 

453.3 

ff 

454.9 

452.5 

frc 

453.2 

450.0 

frcr 

453.1 

449.9 

Data 

453 

1.3 


chamber specific impulses (ISP’s) for 3-D cases 


Table 5 shows the comparison of computed specific impulses for the 3-D cases. The 
qualitative trend among the cases is very similar to the corresponding axisymmetric cases (Table 
4). Hie results of grid 3d9 are consistently lower than those of grid 3d6. This is because the 
effective cell density of grid 3d9 is less than that of 3d6, although the total number of cells in 
grid 3d9 is higher than that of grid 3d6 (Table 1). As a general rule of thumb, the accuracy of 
two tetrahedral cells is approximately equivalent to that of one hexagonal cell. On the other 
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hand, since the number of cells in grid 3d9 is more than those in grid 3d6, it costs more to run 
grid 3d9. This demonstrates that the structured-cell dominated grid 3d6 is more favorable both 
in terms of accuracy and computational efficiency, similar to the findings in the axisymmetric 
cases. This also agrees with the result of Huynh’s Fourier analysis 13 that the upwind scheme 
prefers structured meshes. Within grid 3d6, again the result of case frc compares very well with 
that of the measurement, while the addition of radiation (case frcr) changes the value only 
slightly. 


VI. Conclusions 

Unified computational analyses for computing the design parameters such as the axial thrust, 
convective and radiative wall heat fluxes for hydrogen-fueled liquid rocket engine thrusters were 
conducted, in order to develop a computational strategy for computing those design parameters 
through parametric investigations. The computational methodology is based on a multi- 
dimensional, finite-volume, turbulent, chemically reacting, radiating, unstructured-grid, and 
pressure-based formulation. Systematic parametric studies on effects of wall boundary 
conditions, combustion chemistry, radiation coupling, computational cell shape, and grid 
refinement were performed and assessed. Under the computational framework of this study, it is 
found that the structured-mesh performed more favorably than the unstructured-mesh. The 
effect of radiation coupling was shown to not make an appreciable difference, while that of grid 
refinement sharpens shock capturing. Finite-rate chemistry option performed better than that of 
the equilibrium chemistry, while the frozen chemistry option is undesirable, due to 
thermodynamics considerations. For regeneratively cooled engines, incorporating the effect of 
heat loss drastically improves the axial force predictions. The computed flow physics, axial 
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thrust performance, and wall heat fluxes compared well with those of available test data and 
design calculations, when the desired computational strategy (structured-grid dominated mesh, 
finite-rate chemistry, and cooled wall) was used. 
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Fig. 1 The layout of hybrid grid axl. Top: the overall 
grid. Bottom left: close-up near the throat. Bottom 
right: close-up near the nozzle lip. 



Fig. 3 Layout of hybrid grid 3d6. Upper figure: an 
overall view. Lower left: a cross-sectional cut through 
the axis. Lower right: the exit plane. 






Fig. 2 The layout of hybrid grid ax6. Top: the overall Fig. 4 Layout of hybrid grid 3d9. Lower left insert: 

grid. Bottom left: close-up near the throat. Bottom the cross-sectional cut of the thruster inlet Lower 

right: close-up near the nozzle lip. right insert: the cross-sectional cut of the exit plane. 
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Fig. 5 A comparison of computed and specified 
regeneradvely cooled wall temperatures for grid ax6. 


Fig. 7 A comparison of computed centerline H 2 O mass 
fractions, specific heat ratios and Mach numbers for 


Fig. 6 Computed Mach number contours, a) case fe, 
grid axl; b) case frcgr, grid axl; c) case fz, grid ax6; 
and d) case frcgr, grid ax6. 


Fig. 8 A comparison 
temperatures for grid ax6. 
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